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1.  Introduction 

The  purpose  of  this  report  is  to  review  develppments  in 
objective  analysis  of  meteorological  fields.  Tlis  first  part 
will  review  the  process  known  as  "optimvun  interpolation”  and 
observe  that  it  coincides  with  schemes  developed  Independently 
in  other  scientific  disciplines.  The  term  "optimum"  arises 
from  the  fact  that  the  expected  mean  squared  error  over  some 
ensemble  of  realizations  (e.g.,  over  time)  is  minimized.  The 
term  "interpolation"  is  misleading  since  it  refers  to  inferring 
values  at  points  other  than  data  points,  but  not,  however, 
by  a  scheme  that  necessarily  reproduces  given  values  at  the 
data  points.  Because  of  this  fact,  it  would  probably  be 
better  to  call  the  process  "optimum  approximation".  However, 
we  follow  the  meteorological  literature  and  retain  the  term 
"optimum  interpolation". 

The  second  basic  thrust  of  this  report  is  to  discuss 
Cressman's  scheme  of  successive  approximations  and  show  that 
a  certain  variant  of  the  scheme  will  converge  to  the  same 
result  given  by  optimum  interpolation.  Use  of  this  process 
could  be  advantageous  from  a  computational  viewpoint,  compared 
to  optimum  interpolation. 

In  Section  2  we  derive  the  optimum  interpolation  scheme 
and  show  the  functional  form  of  the  approximation.  We  address 
some  computational  aspects  and  recent  developments  in  Section 
3.  In  Section  4  Cressman's  successive  correction  scheme  is 
discussed,  including  a  statistically  motivated  variant  of  it. 
The  final  section  is  devoted  to  showing  that  a  suitable 


variant  of  a  successive  corrections  scheme  will  converge  to 
the  same  function  as  given  by  optimum  interpolation. 

2.  The  Functional  Form  of  Optimum  Interpolation 

The  develofanent  of  optimm  interpolation  in  meteorology 
dates  back  to  Gandin  [12]  emd  is  based  on  Wiener-Kolmogorov 
theory  in  time  series  analysis.  Various  disciplines  have  used 
similar  schemes  for  some  time,  apparently  developed  indepen¬ 
dently.  We  have  discovered  references  to  developments  in 
geology /mining  (where  it  is  usually  called  kriging)  [1] ,  [20] , 
[22],  [23],  [27],  photogrammetry  [19],  geodesy  [15],  statistics 
(where  it  is  called  stochastic  process  prediction)  [32] »  and 
electrical  engineering  (where  it  is  sometimes  called  a  Wiener 
filter)  [8]. 

We  derive  the  general  form  of  optimum  interpolation  and 
show  the  foimi  of  the  interpolation  function.  While  the  latter 
is  known,  it  is  not  well  known  and  is  even  disavowed  in  print 
in  one  paper  [reply  to  2].  This  situation  has  probably  occurred 
because  the  principal  interest  is  to  obtain  a  grid  of  points 
from  scattered  observations  and  not  to  obtain  an  approximating 
surface.  However,  the  form  of  the  equation  of  the  surface  is 
interesting  and  revealing. 

Let  X  be  the  independent  variable,  and  Z(X)  be  a  random 
fvinction  whose  value  is  to  be  estimated  from  known  or  measured 
values  Z(Xj^),  Z(X2),  ...,  Z(^)  at  scattered  points,  Xi»  •••» 

2^.  We  denote  the  expected  value  of  Z (X) ,  E[Z(X)]  by  m(X) . 

This  mean  value  as  a  function  of  position  is  called  the  trend 
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surface,  and  depending  on  the  source  of  the  data,  nay  be 
assumed  to  have  a  particular  form,  or  to  be  zero.  We  will 
assume  that  the  trend  surface  is  given  by 


n 

m(X)  =  I  c.f.  (X)  ,  (n  <  N-1) 
k»0  ~ 

where  the  known  linearly  independent  functions 

(with  unknown  coefficients,  Cj^)  .  For  meteorological  applica¬ 
tions,  Z  represents  a  residual  (deviation  from  climatology, 
say)  which  is  assximed  to  have  zero  mean,  and  thus  m(X)  =  0. 

We  include  the  term  for  completeness  in  our  development. 

A  number  of  assximptions  will  be  made  concerning  the  dis¬ 
tribution  of  the  random  function  Z (X) .  We  want  to  estimate 
Z(X)  by  a  linear  predictor, 

N 

Z(X)  =  I  X.Z(X.)  . 

j=l  3  -3 

We  ass\ime  that  the  optimum  predictor  is  linear  in  the  observed 
values,  which  is  the  case  if  the  distribution  is  Gaussian,  but 
not  necessarily  otherwise.  In  principle  the  following  process 
can  formally  be  carried  out  without  assximptions  about  the 
covariance  function 

(1)  C(X,Y)  =  El(Z(X)  -m(X))(Z(y)  -  m(Y))]  . 

In  practice,  and  for  computational  reasons  it  is  convenient 
to  make  the  assumptions  of  stationarity  and  isotropy  for  the 


'•  .1 v"^V  ' 


covariance  function.  The  net  effect  of  these  assumptions  is 
that  the  covariance  function  C(X,Y)  is  a  function  of  the  distance 
between  X  and  Y  only,  not  X,Y,  or  their  relative  positions 
other  than  distance  between  them.  These  assumptions  probably 
do  not  hold  in  meteorological  applications;  for  example, 
prevailing  winds  will  certainly  tend  to  give  a  distortion  from 
isotropy  and  various  landforms  will  give  a  distortion  from 
stationarity. 

We  want  to  estimate  the  value  of  Z  at  X;  let  us  call  that 
estimate  Z(X).  We  will  do  this  by  minimizing  E[(Z(X)  -Z(X))^] 
where 


Z(X)  = 


I  A.Z(X.)  , 
j=i  3  : 


STibject  to  some  conditions  which  guarantee  unbiased  estimates. 

For  example  if  Z (X)  has  an  unknown  constant  mean,  E(Z(X))  ^  o  , 

N 

then  the  constraint  ^  X .=  1  is  needed  to  guarantee  the  estimate 

j=l  3 

is  unbiased.  In  the  general  case,  the  constraints  to  be  im¬ 
posed  are 


j  ^  j  K 


k  ~  0,1, ...,n« 


Note  that  this  implies  we  must  have  n  <  N  and  further  it  can 
be  deduced  that  if  the  data  lies  on  the  trend  surface,  so  will 
the  estimated  point  (provided  the  estimate  is  xinique,  a  standard 
assumption) . 
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In  meteorological  applications  it  is  assumed  the  measure¬ 
ments  ("known"  values,  Z (X^) )  are  subject  to  errors,  hence  the 
measured  values  are  Z (X^)  +  £(X^).  We  assume  the  errors  are 
Gaussian  with  mean  zero,  and  are  independent  of  the  function 
Z (X) .  We  denote  the  covariance  function  for  the  errors  by 
C^(X,Y).  Ultimately  we  assijme  the  errors  are  independent, 
so  that  we  will  have  C  (X,Xj)  =  6(X-X.),  where  6(0)  =  1, 

6  (X)  =  0,  X  0.  For  the  derivation,  however,  we  will  allow 
the  more  general  covariance  function.  We  note  that  in  the  case 
of  satellite  data,  for  example,  the  assumption  of  independence 

and  zero  mean  will  probeUaly  not  be  satisfied. 

~  2 

To  minimize  E ( (Z (X)  -  Z (X) )  ]  subject  to  the  constraints 

(2)  ,  we  use  Lagrange  multipliers,  2yj^,  obtaining  the  objective 
function 

N  ,  n  N 

(3)  EI(Z(X)-  I  A.  (Z(X.)  +e(X  )))^]  +  I  2y  (  I  A  f  (X  )  -  f (X)  ) 

j=l  J  k=0  j=l 


Before  differentiation,  we  write  this  as 


(4)  E[(Z(X)  -m(X))^  -  2(Z(X)  -m(X))  A  .  (Z  (X.)  +  e  (X.)  -  m(X. )  ) 

-  -  -  3  3 


N 

+  (  y  A.  (Z(X.)  +  e(X.)  -m(X.)) 
j_l  33  3  3 


Upon  taking  partial  derivatives  with  respect  to  the  and 
yj^,  and  simplifying  somewhat,  we  obtain 


I  >jlc(x.,x.)  +  c^(Xi,Xj)l  +  i‘^fk(5i) 


C(x,x^) 


X  —  1,2, •••,!(, 


and  the  constraint  egxiations  (2) 


In  matrix  form  the  system  of  (linear)  equations  to  be 
solved  is  conveniently  represented  in  partitioned  form  as 


M  ;  F 

I  0 


where 


M  =  (C(X^,Xj)  +  C^CX^rXj))  iO  =  lr..wN 


F  —  i— 


0  is  a  zero  matrix  of  order  (n+1)  x  (n+1) 


_X  -  (Xj^,X2/ .  .  . /Ajj) 


=  (C(X,X^),  ...f  C(X,3^))^,  and 


=  (f-.(X).  f_(X))''. 


Letting  G  denote  the  coefficient  matrix,  we  have  that 
(formally)  the  solution  is 


Letting  d  «  (Z  (X^^)  +  e  (X^^) , . . . , Z (^)  +  e  (3^)  ,0, . . .  ,0)  represent 


the  data  vector,  we  obtain 


Z  (X)  = 


I  X.  (Z(X.)  +  e(X.) 
j=l  ^ 


=  d‘ 


=  d^G-1 


We  will  give  an  alternative  interpretation  of  Z.  Note  that 


depends  on  the  value  of  X  as  well  as  the  data  points 


Xj^,,..,3^,  while  d  is  independent  of  X.  Since  G  is  symmetric, 
so  is  g”^,  and  we  have 


Z(X)  =  d^G"^ 


=  (G”^d)^( - )  . 


Now,  G”^d  is  the  solution  of  a  certain  system  of  equations, 
ncunely 


Ga  —  d. 


where  a  =  bQ...bj^)*'.  This  represents  Z(X)  in  the 
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form 


(9) 


“  ■ 

N  n 

=  I  AiC(X,X.)  +  I  b.f.  (X) 
i=l  ^  ^  k=0  *  ^  “ 


The  system  of  equations  (8)  can  be  thought  of  as  arising  from 
the  requirements  that  an  approximation  consisting  of  a  linear 
combination  of  the  functions  CCX^X^)  +  C^(X,X^),  i  = 
and  f^(X),  k  »  0,...,n  be  required  to  interpolate  the  data, 
Z(X^)  +  e  (Xj^)  ,  i  =  1,...,N,  along  with  constraints  analogous 
to  exactness  for  the  f]j(X), 

N 

\  A . f .  (X • )  —  0,  k  —  0, !,•••, n  • 

j=l  J  “J 

N 

Of  course  the  terms  I  A.C  (X,X. )  represent  interpolation  to 

i=l  ^ 

the  error  function  and  are  then  dropped  to  obtain  (9) .  View¬ 
ing  things  from  this  perspective,  the  computation  of  Z  at  a 
number  of  different  points  is  simplified,  provided  the  error 
estimate  given  by  optimum  interpolation  is  not  to  be  computed. 
We  address  this  briefly  in  the  next  section. 

The  point  of  view  afforded  by  (9)  makes  it  apparent  that 
"regionalizing"  the  process  by  choosing  (from  a  larger  set) 
data  points  near  the  X  of  interest  must  lead  to  a  discontinuous 
surface  which  may,  in  turn,  lead  to  unne  jessary  and  unwanted 
disturbances.  Phillips  [31]  addresses  this  problem  when  dis¬ 
cussing  combined  analysis  and  initialization  (or  perhaps. 


better  to  say,  analysis  which  does  not  require  initializa¬ 
tion)  .  See  also  Williamson  and  Daley  for  an  iterative  approach 
to  overcoming  this  problem. 

In  meteorological  applications  the  error  covariances  are 

assumed  independent  [see,  e.g.,  7],  in  which  case  C^(X,X^)  = 

2  2 

a  6(X-X.),  where  a  is  the  variance  of  the  error  at  X.. 

^i 

In  this  instance  the  matrix  M  differs  from  the  matrix  (C(X,Xj)) 
only  in  that  the  diagonal  terms  are  augmented.  This  has  a 
beneficial  effect  in  terms  of  the  condition  number  of  G  and 
hence  the  numerical  process  of  solving  (6)  or  (8) . 

The  equations  we  have  derived  are  for  a  single  dependent 
variable  Z (X) .  In  both  meteorology  and  geology  simultaneous 
treatment  of  related  dependent  variables  has  been  derived  and 
is  used.  If  all  dependent  variables  are  measured  at  each  data 
point,  X^,  and  if  the  sum  of  the  expected  squared  errors  for 
the  dependent  estimates  is  to  be  minimized,  the  final  result 
is  formally  the  same  as  equation  (6) .  However,  each  entry  in 
M,  and  each  of  Z (X^)  and  Xj  must  be  vectors,  and  the  entries 
in  the  matrix  F  are  replaced  by  block  diagonal  matrices.  See 
Myers  [24]  for  details,  where  the  process  is  called  co-kriging. 

In  meteorology  not  all  variables  are  measured  at  each  data 
point.  The  complication  this  causes  is  readily  resolved, 
although  it  is  simpler  to  group  variables  instead  of  points. 

It  is  called  multivariate  optimum  interpolation  [7],  a  confusing 
term  to  those  outside  the  field  since  the  "multivariate" 
refers  to  the  dependent  variables,  not  the  independent  variables 
Of  course,  cross  covariances  between  the  dependent  variables 


are  required.  See  [7]  and  [21]  for  the  development  In 
meteorological  terms. 

As  a  matter  of  interest,  we  observe  that  the  process  is 
somewhat  reminiscent  of  Lagrange  interpolation,  with  the 


playing  the  part  of  the  fundamental  Lagrange  polynomials. 
Thus,  solution  of  (6)  for  the  Xj  is  equivalent  to  solving  for 
the  values  of  the  fundamental  Lagrange  polynomials  at  the 
point  X.  Alternatively,  solving  (8)  for  the  A^  is  equivalent 
to  solving  for  coefficients  in  the  interpolation  polynomial 
expressed  as  a  linear  combination  of  polynomials. 


3.  Practical  Considerations  and  Recent  Results 


One  of  the  most  important  aspects  of  optimum  interpolation 
is  the  appropriate  specification  of  the  covariance  function, 
C(X/Y).  In  meteorology  this  has  been  treated  by  a  number  of 
authors,  [3],  [9],  [33].  The  importance  of  this  has  been 
recently  noted  by  Franke  [11],  Hollingsworth  [16],  and  Lorenc 
[21]  from  a  practical  point  of  view.  In  theory,  Yakowitz  and 
Szidarovszky  [43]  have  shown  that  (in  the  absence  of  measure¬ 
ment  errors)  the  approximation  converges  as  the  set  of  data 
points  becomes  dense.  Within  some  limitations  this  result 
holds  even  if  the  covariance  functions  are  wrong. 

The  error  estimate  for  optimum  interpolation  can  be  shown 
(by  substitution)  to  be 


E[(Z(>’  -  Z(X)  )“"]  = 


C(X,X)  -  V^(2Q^m"^  -  Q^m"V)V^  , 
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where  Q  =  I  -  F  (F^m“^F)  ”  .  If  E[Z(X)]  =  0,  then  F  =  0 

and  the  above  reduces  to  the  more  familiar  form 

E[(Z(X)  -  Z(X)  )^]  =  C(X,X)  -  . 

As  noted  by  Yakowitz  and  Szidarovsky  [43] ,  this  estimate  is 
good  only  if  the  covariance  functions  are  correct.  They  show 
that  error  estimates  with  incorrect  covariance  functions  may 
be  so  poor  as  to  not  converge  to  zero  as  the  data  points  be¬ 
come  dense,  even  though  the  approximation  converges.  The 
net  result  of  this  is  that  one  should  not  place  too  much  faith 
in  the  error  estimates.  The  covariances  assumed  are  almost 
certainly  wrong,  and  the  more  drastic  effect  is  on  the  error 
estimate  rather  than  the  approximation  itself. 

Computationally  the  choice  between  solving  (6) ,  then 
evaluating  Z (X)  by  (7),  or  solving  (8),  then  evaluating 
Z (X)  by  (9)  depends  on  two  things:  (i)  If  the  error  estimate 
is  also  to  be  computed,  (6)  -  (7)  is  cheaper;  (ii)  If  the 
error  estimate  is  not  to  be  computed,  then  (8)  -  (9)  is  cheaper, 
except  in  the  instance  of  only  one  evaluation  of  Z (X) . 

In  meteorological  applications  it  is  impossible  to 
sider  all  data  points  at  once.  This  leads  to  some  selection 
process  based  on  the  "most  important"  observations.  Often  the 
closest  points  are  considered  the  most  important.  Another  scheme 
is  to  retain  the  points  corresponding  to  the  larger  terms  in  the 
covariance  matrix  (C(X^,Xj)).  Since  the  importance  of  a  point 
depends  on  the  entries  in  the  inverse  of  the  matrix  rather 
than  the  matrix  itself,  this  does  not  seem  to  be  a  good  scheme. 
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A  reasonable  choice  is  probably  to  use  "closest  points"  in 
some  norm  which  accovints  for  prevailing  influences r  e.g., 
winds . 

Several  recent  papers  of  practical  and  theoretical  inter¬ 
est  relate  optimum  interpolation  to  conventional  approximation 
theory.  Among  these  are  Kimmeldorf  and  Wahba  [17],  Matheron 
[24],  Salkauskas  [33],  and  Wahba  [40].  The  significance  of 
the  result  for  meteorological  applications  is  discussed  by 
Wahba.  For  errors  which  have  common  variance  and  covariance 
functions  which  have  a  finite  square  integral  (over  the  entire 
plane) ,  optimum  interpolation  leads  to  the  function  Z (X)  given 
by  (7) .  This  function  is  also  the  solution  of  a  variational 
problem  in  the  spatial  domain,  that  is;  Find  Y(X)  (in  a  cer¬ 
tain  reproducing  kernel  Hilbe.vt  space)  to  minimize 

N  P 

I  (Z(X.)  +  e  -  Y(X.))"^  +  XJ(Y)  , 
j=l  3  D  -3 

where  X  is  the  ratio  of  the  variance  of  the  errors  to  the 
variance  of  the  random  function  Z(X)  (=  C(X,X)),  and  J(Y) 
is  the  square  norm  of  Y  in  the  Hilbert  space.  By  appropriate 
choice  of  the  functional  J(Y) ,  new  methods  can  be  obtained 
which  minimize  or  eliminate  contributions  from  unwanted  modes. 

4.  Cressman's  Successive  Correction  Scheme 


converge  (although  as  applied,  usually  does)  if  the  iteration 
is  continued.  In  the  next  section,  we  will  show  a  relation 
between  a  variant  of  the  Cressman  scheme  and  optimum  inter¬ 
polation.  We  note  that  Cressman 's  scheme  bears  a  resemblance 
to  Shepard's  method  [35],  [13],  but  its  differences  are  more 
important  than  its  similarities.  Most  importantly,  the  gradi¬ 
ents  of  the  approximating  surface  are  not  necessarily  zero  at 
the  data  points  as  they  are  for  Shepard's  method.  In  addition, 
as  originally  proposed  by  Cressman,  the  function  is  not  smooth, 
i.e.,  does  not  have  continuous  partial  derivatives. 

Cressman 's  scheme  achieves  a  weighted  average  of  the  data 
(a  convex  combination,  in  fact)  as  follows.  Let  the  data  again 
be  denoted  by  Z(Xj),  j  =  1,...,N,  and  associate  a  weight  func¬ 
tion,  Wj (X) ,  with  each  point  X^ .  This  function  is  ordinarily 
a  univariate  function  of  distance,  l|X-Xjl|.  Cressman 
proposed  using 


(10) 


Wj  (X) 


+  D 


Ix-x.ll^  < 


0  , 


IX-X.  M  >  D‘ 


Another  scheme  [4]  is  to  take 


(11)  Wj(X)  =  exp(-l  |x-XjM  Vb^)  . 

A  first  approximation  is  taken  to  be 

N  N 

(12)  Z^^MX)  =  I  W.(X)Z(X.)/  1  W.(X). 

j=l  ^  j  =  l  => 


13 


The  denominator  normalizes  the  weights  and  since  the  (X) 

are  positive,  (X)  is  a  convex  combination  of  the  data. 

As  such,  it  satisfies  the  property,  min (Z (X<) )  <  Z^®^  (X)  £ 

j  ^ 

max(Z  (X-)  )  . 
j  ^ 

Usually,  the  process  is  repeated  to  correct  the  differ¬ 
ences  between  the  data  and  the  approximation  Z(Xj)  -  Z^^^  (Xj) , 
but  using  a  smaller  "radius".  This  means  using  a  smaller  D 
in  (10),  or  a  smaller  B  in  (11).  With  superscripts  denoting 
iterations,  we  have  the  following  scheme:  Let  Z^^^(X)be  given 
by  (12),  with  (X)  replaced  by  (X).  Then 

(13)  Z^^^(X)  =  Z^^“^^(X)  +  I  W.  ^^Nx)  (Z(X.) 

j=l  =•  D 

-  Z^^"^Nx.)/  f  wJ^^(X),  k  =  1,2...  . 

“3  -4=1  3  “ 


If  we  look  at  the  sequence  of  vectors  which  approximate  the 
data  vector,  Z=  (Z  (X^^)  . .  .Z  (3^) )  ^,  we  have  {(Z^^Mx^)... 

Z^*^^  (3^))^},  k  =  0,1,...  .  Denote  these  vectors  by  ,  and 

define  the  matrix 


(14) 


h(W 


r«  1 


Then  the  iteration  takes  the  form 


i,j  =  If 


•  i 


N. 


Z  Z 


(15) 


+  H^^Nz  -  )  ,  k  =  1,2,... 


14 


=  (I (Z  -  ,  k  =  1,2,... 

This  easily  leads  to 

(16)  Z-Z^^^  =  [n  (I  -  )  ]  (Z  -  Z^°^  )  , 

P=1  -  - 

and  we  see  that  this  iteration  converges  provided  that,  for 

sufficiently  large  k,  the  norms  of  the  I-H'  '  are  bounded  by 

a  constant  less  than  one.  This  holds  for  any  norm;  hence,  if 

(k) 

all  eigenvalues  of  the  I-H  have  magnitude  bounded  by  a 
constant  less  than  one,  for  sufficiently  large  k,  convergence 
is  obtained. 

Generally,  the  effect  of  decreasing  D  in  (10)  or  B  in  (11) 

is  to  increase  the  relative  size  of  the  diagonal  elements  in 
(k) 

H  .  Since  each  row  sums  to  one,  the  matrix  will  eventually 

become  diagonally  dominant.  In  any  case,  if  all  eigenvalues 

are  positive  (as  for  (11)  ,  for  example) ,  the  eigenvalues  of 
(k) 

I-H  are  then  bounded  by  a  constant  less  than  one,  inde¬ 
pendent  of  k,  provided  B  is  a  decreasing  function  of  k. 

The  situation  for  weights  given  by  (10)  is  not  so  pleasant 

(k) 

In  this  case,  the  matrix  H  may  have  negative  eigenvalues, 

( k) 

which  leads  to  I-H  having  eigenvalues  greater  than  one. 


As  D  decreases,  all  eigenvalues  do  become  positive.  Whether 
or  not  this  happens  for  values  of  D  used  in  practice  should 
be  investigated  since  this  has  a  bearing  on  the  stability  of 
the  iteration. 

The  effect  of  decreasing  D  in  (10)  or  B  in  (11)  is  to  speed 

convergence  of  the  iteration,  since  this  tends  to  increase  the 

(k)  (k) 

eigenvalues  of  H  .  We  note  in  passing  that  H'  '  is  a  sto¬ 
chastic  matrix,  and  thus  has  its  largest  eigenvalue  equal  to 

T 

one,  with  eigenvector  (l,...l)  .  In  terms  of  the  approxima¬ 
tion,  this  means  convergence  in  one  iteration  if  ^  is  a  con¬ 
stant  vector,  i.e.,  (c,...,c)^,  c  is  a  constant.  The 

maximum/minimum  principle  cited  earlier  implies  the  scheme  is 
exact  for  constants  however. 

In  the  next  section  we  discuss  the  general  form  of  the 
approximation,  and  show  that  under  certain  simple  modifications 
the  scheme  approximates  optimum  interpolation. 

5 .  Relation  of  a  Variant  of  Pressman's  Scheme  to  Optimum 
Interpolation 

The  general  form  of  Z^®^(X)  in  (12)  is  a  rational  function 

in  the  weights,  (X) ,  and  if  the  scheme  is  iterated  as  in  (13) 

the  form  is  that  of  a  sum  of  functions,  each  rational  in  the 

(k) 

appropriate  set  of  weights  ,  k=0,l,...  .  If  weights 

were  taken  to  be  the  functions  C(X,Xj)  +C^(X,Xj)  in  (1),  for 
all  k,  the  resulting  approximation  bears  some  relationship  to 
optimum  interpolation,  although  it  is  rational  in  the  covari¬ 


ance  functions  rather  than  linear  in  them.  However,  if  the 

(k) 

denominators  of  (12)  and  (13)  (and  hence,  of  H  )  are  replaced 


by  a  suitable  constant,  the  iteration  will  converge  to  the 
optimum  interpolation  function. 

First  we  consider  the  covariance  matrix 


(17)  M  =  (C(X^,Xj)  +  C^(Xj^,Xj))  ,  i,j  =  1,...,N  . 

This  matrix  must  be  positive  semidefinite  and  we  make  the 
usual  assumption  that  it  is  definite,  i.e.,  has  no  zero  eigen¬ 
values.  Let  I |m| I  denote  the  max  row  sum  norm  of  M,  and  let 
B  be  a  constant  satisfying  B  >  jM] | .  Now  consider  the 
iteration  obtained  by  replacing  the  denominators  in  (12)  and 
(13)  by  B/  which  leads  to  the  matrix  iteration  analogous  to 


2(k)  ^  +  |m(Z  -  ) ,  k=l,2,...  . 


This  leads  to  the  analogue  of  (16) , 


Z  -  Z 


(I  -|m)^(Z  -  Z^°^ 


Thus,  convergence  is  obtained  whenever  has  all  eigen- 

p 

values  strictly  bounded  by  one.  Since  the  eigenvalues  of 

must  be  bounded  by  |  |^M  |  |  <  2  by  our  choice  of  B»  convergence 

(k) 

is  obtained.  The  form  of  each  Z'  MX)  is  a  linear  combination 

of  the  C(X,X.)  +  C  (X,X.) .  Because  convergence  implies  agree- 
J  e  —  3 

ment  at  the  data  points,  we  see  that  if  the  error  covariances 
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have  the  form  noted  before,  C  (X,X. )  »  6(X-X.),  the 

C  ^  "“i 

limiting  approximation  given  by  (18)  agrees  exactly  with  that 
given  by  (8)  -  (9) ,  except  at  the  X^,  where  a  jump  occurs  to 
yield  "interpolation."  Of  course  one  thinks  in  terms  of 
dropping  that  term  for  the  final  approximation,  but  must  do 
so  only  if  an  evaluation  point  coincides  with  a  grid  point. 
Practically  speaking,  the  reverse  situation  is  where  the 
special  instance  is  encountered. 

The  rate  of  convergence  may  be  slow  because  of  the  like¬ 
lihood  of  M  having  small  eigenvalues,  leading  to  I-^M  having 

P 

eigenvalues  close  to  one.  However,  the  presence  of  the  error 
covariance  tends  to  increase  the  eigenvalues  of  M,  and  in  this 
respect,  large  observational  errors  would  benefit  the  conver¬ 
gence  rate.  It  would  seem  best  to  try  to  choose  3  to  minimize 
the  magnitude  of  the  eigenvalue  of  I-^M  of  largest  magnitude. 
This  would  maximize  the  rate  of  convergence.  On  the  other 
hand,  most  of  the  significant  information  may  correspond  to 

large  eigenvalues  of  M.  (Recall  that  one  is  an  eigenvalue  of 

(k)  t 

H  With  eigenvector  (1,1,.,.,1)  .)  In  this  case,  it  would 

make  sense  to  take  3  ~  ||m||  which  would  cause  rapid  conver¬ 
gence  for  these  modes,  while  modes  corresponding  to  small 
eigenvalues  are  of  high  frequency  and  could  be  best  filtered 
out. 

The  filtering  potential  of  this  scheme  should  be  investi¬ 
gated  further  to  determine  whether  or  not  the  eigenvectors 
corresponding  to  small  eigenvalues  do  indeed  lead  to  unwanted 
noise  in  the  approximation  which  later  must  be  filtered  out. 
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If  so,  this  scheme  could  be  an  advantageous  one  to  use. 

Some  simulations  of  the  scheme  have  been  carried  out  through 
the  iteration  process.  However,  the  results  are  nontrivial 
to  interpret  and  need  additional  study,  particularly  in  the 
light  of  the  filtering  scheme  presently  used  in  the  opera¬ 
tional  model.  The  combination  of  including  the  measurement 
errors  and  the  constant  normalization  factor  will  result  in 
the  successive  correction  method  appearing  more  like  optimum 
interpolation.  A  multivariate  scheme  could  be  derived  in 
a  straightforward  fashion. 
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31.  Director 

Geophys .  Fluid  Dynamics  Lab 
NOAA,  Princeton  University 
P.O.  Box  308 
Princeton,  NJ  08540 

32.  Head,  Atmos.  Sciences  Diy. 

National  Science  Foundation 
1800  G  Street,  NW 
Washington,  DC  20550 

33.  Laboratory  for  Atmos.  Sciences 
NASA  Goddard  Space  Flight  Cen. 
Greenbelt,  MD  20771 

34.  Atmospheric  Sciences  Dept. 

UCLA 

405  Hilgard  Ave. 

Los  Angeles,  CA  90024 
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Director 

Institute  of  Geophysics 
Univ.  of  Calif.  Los  Angeles 
Los  Angeles,  CA  90024 

Chairman 

Meteorology  Dept. 

California  State  University 
San  Jose,  CA  95192 

National  Center  for  Atmos. 

Rsch. ,  Library  Acquisitions 
P.O.  Box  3000 
Boulder,  CO  80302 

Director 

Denver  Research  Institute 
University  of  Denver 
3115  S.  University  Blvd. 
Denver,  CO  80210 

Colorado  State  University 
Atmospheric  Sciences  Dept. 
Attn;  Dr.  Roger  Pielke 
Fort  Collins,  CO  80523 

Nova  University 
Physical  Oceanographic  Lab 
8000  North  Ocean  Dr. 

Dania,  FL  33004 

Chairman 

Meteorology  &  Physics  Dept. 
University  of  Florida 
215  Physics  Bldg. 

Gainesville,  FL  32601 

University  of  Miami 
R.S.M.A.S.  Library 
4600  Rickenbacker  Causeway 
Virginia  Kay 
Miami,  FL  33149 

Florida  State  University 
Environmental  Sciences  Dept. 
Tallahassee,  FL  32306 

University  of  Hawaii 
Meteorology  Dept. 

2525  Correa  Road 
Honolulu,  HI  96822 


45.  Geophysical  Sciences  Collection 
The  University  of  Chicago 

The  Joseph  Regenstein  Library 
1100  East  57th  Street 
Chicago,  IL  60637 

46.  University  of  Maryland 
Meteorology  Dept. 

College  Park,  MD  20742 

47.  Johns  Hopkins  University 
Applied  Physics  Laboratory 
R.E.  Gibson  Library 
Johns  Hopkins  Road 
Laurel,  MD  20810 

48.  Chairman 
Meteorology  Dept. 

Massachusetts  Institute  of 

Technology 

Ccunbridge,  MA  02139 

49.  Chairman 

Meteorology  &  Oceano.  Dept. 
University  of  Michigan 
4072  E.  Engineering  Bldg. 

Ann  Arbor,  MI  48104 

50.  Chairman 
Physics  Dept. 

University  of  Minnesota 
Minneapolis,  MN  55455 

51.  Director  of  Meteorology 
Earth  &  Atmos.  Science  Dept. 

St.  Louis  University 

P.O.  Box  8099 

St.  Louis,  MO  63156 

52.  Atmospheric  Sciences  Center 
Desert  Research  Institute 
P.O.  Box  60220 

Reno,  NV  89506 

53.  Atmospheric  Sci.  Rsch.  Center 
New  York  State  Univ. 

1400  Washington  Ave. 

Albany,  NY  12222 

54.  Reference  Library 
Maritime  College 
New  York  State  Univ. 

Fort  Schuyler 
Bronx,  NY  10465 
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New  York  Polytechnic  Institute 
Meteoro.  &  Oceano.  Dept. 

333  Jay  Street 
Brooklyn,  NY  11201 

Chairman 

Meteorol.  &  Physical  Oceano. 
Cook  College,  P.O.  Box  231 
Rutgers  University 
New  Brunswick,  NJ  08903 

Chairman 

Meteorology  Dept. 

University  of  Oklahoma 
Norman,  OK  73069 

Professor  Y.  Sasaki 
University  of  Oklahoma 
Rm.  219,  200  Felgar  St. 

Norman,  OK  73069 

Atmospheric  Sciences  Dept. 
Oregon  State  University 
Corvallis,  OR  97331 

Dean  of  The  College  of  Science 
Drexel  Institute  of  Technology 
Philadelphia,  PA  19104 

Chairman,  Meteorology  Dept. 
Pennsylvania  State  Univ. 

503  Deike  Bldg. 

University  Park,  PA  16802 

Texas  A&M  University 
Meteorology  Dept. 

College  Station,  TX  77843 

Director  of  Research 
Institute  for  Storm  Research 
University  of  St.  Thomas 
3812  Montrose  Blvd. 

Houston,  TX  77006 

Chairman 

Meteorology  Dept. 

University  of  Utah 

Salt  Lake  City,  UT  84112 

Chairman 

Atmos.  Sciences  Dept. 

University  of  Virginia 
Charlottesville,  VA  22903 
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.  Director 

Oceanographic  Institute 
Old  Dominion  University 
Norfolk,  VA  23508 

67.  University  of  Washington 
Atmospheric  Sciences  Dept. 
Seattle,  WA  98195 

68.  Chairman,  Meteorology  Dept. 
University  of  Wisconsin 
Meteoro.  &  Space  Science  Bldg. 
1225  W.  Dayton  St. 

Madison,  WI  53706 

69.  Colorado  State  University 
Atmospheric  Sciences  Dept. 
Attn;  Librarian 

Fort  Collins,  CO  80521 

70.  Saint  Louis  University 
Earth  &  Atmos.  Sciences  Dept. 
P.O.  Box  8099-Laclede  Station 
St.  Louis,  MO  63156 

71.  Woods  Hole  Oceanographic  Inst. 
Document  Library  LO-206 
Woods  Hole,  MA  02543 

72.  Colorado  State  University 
Atmospheric  Sciences  Dept. 
Attn;  Dr.  William  Gray 
Fort  Collins,  CO  80523 

73.  Control  Data  Corp. 

Meteorology  Dept.  Rsch.  Div. 
2800  E.  Old  Shakopee  Rd. 

Box  1249 

Minneapolis,  MN  55440 

74.  Systems  &  Applied  Sci.  Corp. 
Attn;  Library,  Suite  500 
6811  Kenilworth  Ave. 

Riverdale,  MD  20840 

75.  Science  Applications,  Inc. 

415  Dela  Vina  Ave. 

Monterey,  CA  93940 

.  Library 

Csiro  Div.  Atmospheric  Physics 
Station  Street 
Aspendale,  3195 
Victoria,  Australia 
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77.  Librarian,  Meteorology  Dept. 
University  of  Melbourne 
Parkville,  Victoria  3052 
Australia 

78.  Library,  Australian  Ntimerical 

Meteorology  Research  Center 
P.O.  Box  5089A 
Melbourne,  Victoria,  3001 
Australia 

79.  Chairman,  Meteorology  Dept. 
McGill  University 

805  Sherbrooke  St.,  W. 
Montreal,  Quebec 
Canada  H3A  2K6 

80.  Meteorological  Office  Library 
London  Road 

Bracknell,  Berkshire 
RG  12  ISZ,  England 

81.  Department  of  Meteorology 
University  of  Reading 

2  Earlygate,  Whiteknights 
Reading  RG6  2AU 
England 

82.  European  Centre  For  Medium 

Range  Weather  Forecasts 
Shinfleld  Park,  Reading 
Berkshire  RG29Ax,  England 

83.  Meteorologie  Nationale 
1  Quai  Branly 

75,  Paris  (7) 

France 

84.  Meteorology  Department 
Andhra  University 
Waltair,  India  530-003 

85.  Hebrew  Univ.  of  Jerusalem 
Library,  Atmos.  Sciences  Dept. 
Jerusalem,  Israel  91904 

86.  Library 

Meteorological  Rsch.  Institute 
1- 1 ,  Nagamine ,  Ya tabe-Machi , 
Ts\ikuba-Gun 
Ibaraki-Ken  305  Japan 
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